Abstract In this article we present a model of formation of a galaxy with 
a black hole in the center. It is based on the Lemaitre-Tolman solution and 
is a refinement of an earlier model. The most important improvement is the 
choice of the interior geometry of the black hole allowing for the formation of 
Gyrs old black holes. Other refinements are the use of an arbitrary Friedmann 
model as the background (unperturbed) initial state and the adaptation of 
the model to an arbitrary density profile of the galaxy. Our main interest was 
the M87 galaxy (NGC 4486), which hosts a supermassive black hole of mass 
3.2-10 9 M Q . It is shown that for this particular galaxy, within the framework of 
our model and for the initial state being a perturbation of the vlCDM model, 
the age of the black hole can be up to 12.7 Gyrs. The dependence of the model 
on the chosen parameters at the time of last scattering was also studied. The 
maximal age of the black hole as a function of the fl m and Qa parameters for 
the M87 galaxy can be 3.717 or 12.708 Gyr. 
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1 Introduction 

Recent years brought increasing evidence that large galaxies host massive black 
holes at their centers. The best evidence was obtained for the Milky Way, due 
to the proximity of the galactic center located at the distance Ro = 8.33 ±0.35 
kpc. Analysis of orbits of the S-stars cluster near the galactic core, especially 
the S2 star, allowed to conclude that a mass of 4.31±0.36-10 6 M Q is responsible 
for the apparent movement of the stars [1,2,3,4]. Such mass enclosed in a small 
volume, of radius less than 0.01 pc, can only be a black hole. The existence of 
black holes at the centers of galaxies calls for construction of models describing 
their creation and growth. This work presents improvements to one of such 
models [5]. 

It is thought that the currently existing structures in the Universe (like 
galaxies, clusters of galaxies and voids) evolved out of small inhomogeneities 
that are observed as directional variations (AT '/T « 10~ 5 ) of the temperature 
of the cosmic microwave background (CMB) radiation. The generation of these 
inhomogeneities is a problem only when one insists on using models of the 
Universe that are born spatially homogeneous. In an inhomogeneous model, 
such as the Lemaitre-Tolman (LT) model used here, the Universe emerges 
from the Big Bang being already inhomogeneous. The inhomogeneities are 
generated within the Big Bang and their amplitudes are arbitrary parameters 
that can be adapted to observational constraints. (But it must be stressed that 
the LT model, in which the matter is assumed to be dust, cannot be applied 
to epochs earlier than last scattering. A still more general inhomogeneous 
model that includes pressure gradients must be used for the pre-last-scattering 
epoch.) It has already been proven in an earlier paper [5], using an LT model, 
that with a suitable choice of the velocity profile at last scattering, an object 
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can be created that has the same density profile at present as the M87 galaxy, 
and contains a black hole around the center that has the same mass as the 
black hole in M87. (For an exposition of the method used in Ref. [5] see Refs. 
[6,7] and [8].) The black hole is created because matter around the center 
expands more slowly than farther away, what leads, at a certain moment, to 
the creation of the apparent horizon surrounding a trapped region, and soon 
after to the creation of a Big Crunch singularity at the center. 

In Ref. [5] two particular configurations of the LT model were considered. 
In one, the black hole formed around a pre-existing wormhole within a fraction 
of a second after the Big Bang. In the other, there was no wormhole or black 
hole initially, and it formed during the evolution of the proto-galaxy. The first 
configuration is somewhat exotic and we will not pursue it here. The problem 
with the second configuration was that the implied age of the black hole was 
only 4 x 10 8 years, what is inconsistent with astrophysical implications [9,10, 
11]. The reason for the implied definite age was that the LT model used there 
contained too few free parameters. 

This paper is an improvement over Ref. [5]. The method of constructing 
the model is the same as before. We define a certain density or velocity profile 
at last scattering, taking care that the amplitude remains within the limits 
implied by the measurements of the temperature distribution of the CMB 
radiation. We define another density profile at the present time that agrees 
with the observationally determined density profile of the M87 galaxy [12] and 
contains a black hole at the center, of mass 3.2 ± 0.9 • 1O 9 M , equal to the 
mass of the black hole observed in M87 [13]. The improvement consists in the 
LT model having more free parameters. Thanks to this, the age of the black 
hole is no longer determined and can be up to ~ 12.7 Gyr, in agreement with 
what astrophysics tells us [9, 10, ll]. 1 

Our model has certain weaknesses, of which we are aware, but, just as was 
stated in Ref. [5], we treat it as an exploratory step into a new territory: we 
intend to test a new method that we hope will be improved in the future. The 
deficiencies to be removed in the future are the following: 

1. No real galaxy is spherically symmetric. The model can be said to apply 
approximately to elliptic galaxies (M87 is one), but spherical symmetry 
leads to the next problem listed below. 

2. Real galaxies rotate. In a spherically symmetric model rotation is necessar- 
ily zero. The presence of rotation influences the time scale of evolution, for 
example, by slowing down the collapse it may significantly delay the forma- 
tion of the black hole. However, no exact solutions of Einstein's equations 
are known that would describe matter that is expanding and rotating at 
the same time. All known expanding solutions have zero rotation, all known 
rotating solutions are either stationary or unrealistic for other reasons (for 
more on this see Ref. [8]). So if one wishes to have an evolving configura- 

1 We do not wish to enter any discussion of the method by which the age of the black 
hole in M87 is inferred. We only wish to demonstrate that, whatever that inference is, our 
model can be made consistent with it. 
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tion that can be described by the exact formalism of general relativity, not 
much is left beyond spherically symmetric models. 
3. As was stated in Ref. [5], the perturbation that would evolve into a single 
galaxy would have the diameter of approx. 0.004° in the CMB sky. The 
current best angular resolution of measurements of the temperature fluc- 
tuations of the CMB radiation is 0.2° [14,15]. Consequently, the amplitude 
of fluctuation at this large angular scale does not give us information that 
we need to constrain our model. Lacking any better possibility, just as in 
the earlier paper [5], we took care that the amplitudes of our initial density 
and velocity profiles do not exceed the limits set at 0.2°. 

The remarks above show that our toy model cannot be literally taken as the 
actual model of an existing galaxy. However, it avoids a few other deficiencies 
that could be contemplated: 

1. Once the black hole is formed in an LT model, where pressure is zero and 
all motions are radial, it keeps accreting matter until it swallows up the 
whole mass contained in the region where E < 0. This happens in a finite 
time, which is arbitrarily long in the neighbourhood of E = 0. However, if 
the function E(M) has such a profile that at a certain M = M g E(M g ) = 
and E(M) > for M > M g , then the mass from the region M > M g is 
not accreted onto the galaxy. Thus, M g , which is independent of time, 
can be interpreted as the mass of the galaxy. This is outside the region 
considered in our paper. The geometrical radius of the M = M g surface 
will be expanding as dictated by the E = evolution equation of the LT 
model. For the surface of a galaxy 2 of radius 32 kpc w 10 21 m and mass 
2.4 x 10 12 M Q = 4.8 x 10 42 kg, if it were to evolve by the E = Lemaitre - 
Tolman equation R, t 2 = 2M/R, the current velocity of expansion would be 
approx. 800 km/s w 8.2 x 10~ 4 pc/yr, and constantly decreasing because 
A = 0. This is negligible compared to the error in determining the edge of 
a galaxy. 

2. Each real galaxy is surrounded by vacuum. The cosmological model takes 
over at a considerably larger scale. If one wants to illustrate this situation 
in a toy model like ours, it is enough to match our LT galaxy model to the 
Schwarzschild spacetime at a certain mass M = M g , and then to match 
the Schwarzschild region on the outside to another LT or Friedmann region 
modelling the Universe. Both matchings would take place outside the region 
we consider and would have no influence on what happens inside the galaxy. 
An explicit example of such matchings is given in Ref. [17]. 

The structure of this paper is as follows. Sec. 2 presents the basic properties 
of the Lemaitre-Tolman model. Sec. 3 briefly outlines the method introduced 
in Refs. [6, 7] and emphasises the key elements of the model. In Sec. 4 we de- 
scribe the improvements to the basic model: the more general forms of the free 

2 Data for the M87 galaxy taken from [16], values of the constants c and G from 
http://physics.nist.gov/cuu/Constants/index.html, and the relations between distance units 
from http://www.asknumbers.com/LengthConversion.aspx. All values rounded off. 
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functions in the LT model and the more general (and more realistic) density 
profile of the galaxy at the present time. This section also gives the necessary 
equations for the use of an arbitrary FLRW model as the background at the 
initial instant (in Ref. [5] the background was assumed to be spatially flat). 
Sec. 5 presents the results of application of the model to the M87 galaxy. Evo- 
lution from spatially homogeneous initial profiles of velocity and density to a 
galaxy with a black hole is described. A spatially homogeneous (flat) density 
profile is obviously within the observational limits on density perturbations, 
but may lead to the following problem. After the LT model is uniquely deter- 
mined by the initial and final density profiles, the initial velocity profile that 
caused the condensation can be calculated from the model and may turn out 
to have a too large amplitude. The same may happen with the initial veloc- 
ity profile being flat - the calculated amplitude of the initial density profile 
may turn out to be too large. Graphs shown in Sec. 5 prove that this did not 
happen; all implied profiles are consistent with the observational constraints. 
Sec. 5 also shows how the age of the black hole and the arbitrary functions in 
the LT model change with the initial FLRW model. The last section contains 
summary and conclusions. 



2 The Lemaitre Tolman model - basic properties 

The Lemaitre-Tolman model is a spherically symmetric nonstatic solution of 
the Einstein equations with a dust source [18,19,20]. Its metric is 

= _ l + 2E(r) dr2 R2{t > r){dd2 + Sin2 M ^ 2) ' (1) 

where E(r) is an arbitrary function arising as an integration constant from 
the Einstein equations, R(t, r) is a function satisfying 

R lt 2 = 2E + 2^ + i/LR 2 , (2) 

where A is the cosmological constant and M(r) is another arbitrary function. 
Equation (2), called the evolution equation, is a first integral of the Einstein 
equations, i?(i, r) is called the areal radius and R, t is referred to as the velocity. 
The matter density p(t, r) is 

8ttG 2M 

R?R yr 

When A = 0, equation (2) has three families of solutions depending on the 
sign of E(r). They are as follows: 
E < - elliptic evolution 

M 

R(t,r) = — ^(1 - cost?), 

(4) 

(_2£)3/2 V ' 

V ~ S[Uf] = M ~ 



4 P=-EZE-- ( 3 ) 
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E = - parabolic evolution 



(J-^/ '/;:-■-': . (-1) 
E > - hyperbolic evolution 

M 

( 2E f/2 ( 6 ) 
smh?7 - t] = — — — (t - t B (r)), 

where i] is a parameter and t B (r) is one more arbitrary integration function 
called the bang time. The formulas given above are covariant under arbitrary 
coordinate transformations of the form r' = g(r), which allows for r to be 
chosen freely, what in turn means that one of the three functions E(r), t B (r) 
and M(r) can be fixed at our convenience by the appropriate choice of g. 

From the form of the solution in the elliptic case, Eqs. (4), we can deduce 
that R(t, r) reaches a maximum for each r and then decreases to zero. The 
values of i] for which i?(t, r) vanishes are -q — 0, corresponding to the Big Bang, 
and rj = 2ir corresponding to the Big Crunch at t = tc 

tc(r)=t B (r)+ { _ 2E)3/2 - (7) 

The function tc(r) is called the crunch time function. 

The Friedman-Lemaitre-Robertson- Walker (FLRW) models arise from the 
LT models in the limit 

t B = const, E = -^kr 2 , M = M r 3 , (8) 

where k and M are constants. 



1/3 



2.1 Shell crossings 

It is important to avoid shell crossings, that is loci of R, r = and M, r ^ 0, 
where a constant r shell collides with its neighbour (as a consequence of g rr = 
0). Shell crossings are curvature singularities and cause the density to diverge 
(p — > oo) and change sign to become negative. Such an undesirable behavior 
can be avoided if the shapes of the three arbitrary functions defining the LT 
model are chosen appropriately. The general conditions for avoidance of shell 
crossings are given in [21,22]. 



Formation of Gyr old BH in the centers of galaxies within LT model 



7 



2.2 M(r) as radial coordinate 

It is convenient to use M(r) as the radial coordinate, that is r' = M(r). It 
is possible because M(r) will be a strictly growing function, M, r > 0, in the 
whole region under consideration. Then we have R = R(t,M), p = p(t,M) 
and Eq. (3) becomes 

8ttG 2 6 



3 f M d£ 



from which we find 

where p is the density profile of the galaxy as a function of M, Ri and 
are constants to be determined later (see Sec. 4.3). Hereafter we use M as the 
radial coordinate. 



3 A galaxy with a central black hole in the LT model 

An LT model is generally specified by defining the shapes of the arbitrary 
functions E(M) and is(M), what then allows for the reconstruction of the 
evolution of the model by using Eqs. (4)-(6). However, it is possible to follow 
a different path, that is to calculate these functions for a specific density or 
velocity profiles at any two time instants, as described in [6] and [7]. Moreover, 
there are no limitations on those profiles (except the obvious one of spherical 
symmetry). The class of the LT evolution applying to each M value can gen- 
erally be different, however for the creation of a black hole an LT model must 
necessarily be recollapsing in the region near the center of symmetry M = 0. 

The main idea of the model constructed in [5] is that a spherically symmet- 
ric galaxy with a central black hole can be created from accretion of matter 
onto an initial fluctuation at recombination. This is the most probable scenario 
of structure formation in our Universe [23,24,25,26,27,28]. Accretion in the 
model is achieved by a slower expansion near the center, and a faster one fur- 
ther away. We will now briefly describe the main properties of the previously 
developed model. 



3.1 The basic framework of the model 

The basic framework of the model was discussed in [5] . The model allows for 
the formation of a galaxy with a central black hole by accretion of background 
matter onto an initial fluctuation of density or velocity distribution at the time 
of last scattering of the CMB. The background is assumed to be Friedmannian 
(see Sec. 4.3). The input data provided for the model are the profiles of the 
velocity or the density at last scattering and the current density distribution in 
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the galaxy. In [5] the model was used to create a black hole from a flat (Fried- 
mannian) velocity distribution at last scattering. In this case, the fluctuation 
was due to the corresponding non-flat density profile. The existence of a black 
hole at the center of a galaxy is assumed on observational grounds. The mass 
of the black hole is an input parameter and likewise must be taken from ob- 
servations. The observational density profile of the galaxy was approximated 
by a simple function, which allowed for exact calculations. 

The construction of the model consisted of the following steps: 

1. determining the functions E(M) and £b(M) in the part outside the black 
hole, that is for the mass range M > Mbh, where Mbh is the observation- 
ally determined mass of the black hole in the galaxy, 

2. smooth joining of the functions E{M) and ts(M) in the parts inside and 
outside the black hole i.e. at the boundary M = Mbh- The part inside the 
black hole is arbitrary and was assumed to be a simple E < recollapsing 
model defined by a choice of the bang and crunch time functions, 

3. evolution reconstruction - using the full E(M) and ts{M) functions the 
evolution of R(t, M) and pit, M) was reconstructed using Eqs. (4)-(6). 

In the model the time of the black hole formation is the value of the future 
apparent horizon at its minimum i.e. (iAH+)min- However, in this model the 
crunch time function and the future apparent horizon almost coincide, the 
instant of the creation of the black hole is very close to tc(0) (see Fig. 4 and 
[5] for details). The first application of the model was to the M87 galaxy which 
hosts a supermassive black hole of mass 3.2 • 10 9 M Q [13], with density profile 
being an approximation of Fabricant et. al. [12]. In the first parametrization of 
the black hole's interior geometry the value of (iAH+)min was high, leading to 
a black hole of the age of a few hundred million years, which is in disagreement 
with other estimations of supermassive black hole's age based on the ylCDM 
model [9,10,11]. 

4 Development of the model 

The model allows for modifications which can lead to more realistic predictions. 
Particular attention was paid to arrive at such bang and crunch time functions 
that would give the black hole's age of the order of Gyr. 

4.1 The black hole interior 

The first modification concentrates on the functions is(M) and tc(M) that 
define the black hole's interior, for which, for fundamental reasons, no observa- 
tional data exist. Therefore, we are free to choose any parametrization and the 
only condition is to join it smoothly to the galaxy model. In the most probable 
scenarios of black hole creation, it is formed by a collapse of a massive object. 
We used exponential and hyperbolic functions to parameterize the interior 
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functions t#(M) and tc(M), see Table 1. The chosen forms of the functions 
have 8 parameters: b\, b 2 , 63, b^ for ts(M) and ci, c 2 , c 3l C4 for tc{M). 

A smooth match of the interior and exterior requires the continuity of 
the functions E(M) and is(M) and their derivatives at the boundary M = 
Mbh (the Darmois/Lichnerowicz junction conditions), this is also a sufficient 
condition [29]. Therefore, the following system of equations has to be solved 
for the parameters bi of i_e(M) and Cj of ic(M), i = 1, . . . , 4: 

tB(M B n, b\, b 2 , b 3 , 64) = [£b] m=Mbh , 



( A/jjji . />, - b 2 . !>:■. h 





~dt B ~ 


)- 


dM 



M=M BH 



E(M B H,bi,b 2 ,b 3 ,b 4 , 

Ci,C2,C 3 ,C 4 ) = [^] m=Mbh 

dE 1 / 

— (^M BH , 61,62, b 3 ,b 4 , 



(11) 





"d£" 


ci,c 2 ,c 3 ,c 4 ^ = 


dM 



M=M BH 

where the expressions on the left hand side are the analytical functions de- 
pending on the chosen parametrization and the expressions on the right hand 
side are the numerically found values of is(M) and E{M) for M = Mbh, us- 
ing the method described in [5]. The analytical expression for E(M) is based 
on Eq. (7). However, the assumed bang and crunch time functions have eight 
unknown parameters, whereas from (11) only four parameters can be deter- 
mined. Therefore, the values of the other four parameters have to be supplied. 
To accomplish that we choose any two parameters for ts(M) and any two pa- 
rameters for tc(M), and set their values using random numbers. The difficulty 
here (from the numerical point of view) is to determine the order of magni- 
tude of the parameters, but this can be achieved by examination of the type 
of parameterized equation (see column 3 of Table 1). The system of equations 
(11) has to be solved numerically. The usage of random numbers does not have 
any impact on the analysis as the multidimensional minimalization methods, 
employed further for finding extreme LT models take the random values only 
as the starting values for subsequent optimalization. 



Name 


t{M) 


Pi 


exp 


Pi + P2 exp (p-iM + pi) 


Pi ~ 

P2 


t(M Bli ) 
~ 10° 


hyp 


Pi + P2 cosh (pzM + pi) 


P3 
Pi 


~ 10 2 
~ 10 2 



Table 1 Parametrizations of the interior bang and crunch time functions. The table shows 
the bang and crunch time functions inside the black hole. The third column gives a general 
order of magnitude of the value of each parameter. t(M^n) means that the value of the 
parameter is of the order of the value of the bang or crunch time function at M = Mbh. 
For ts{M) the parameters pi are denoted bi, and for tc(M) by c^. 
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4.2 The density profile of the galaxy 

A galaxy in the constructed model is described by a density profile and the 
mass of the black hole. However, the quantity used in the numerical calcula- 
tions is the areal radius R(t ,M) corresponding to that (spherically symmet- 
ric) density profile, where to is the present age of the Universe. The equation 
used for finding its value is (10), with p being a function of M . Therefore, the 
general r coordinate representing the distance from the center of symmetry 
must be replaced by the mass M of the sphere of radius r. In [5] the relation 
r = r{M) was found analytically, using an analytically convenient approxima- 
tion to the observational density profile of the galaxy. However, we would like 
to apply the model to any galaxy with known density profile which is suspected 
of hosting a black hole of known mass. In such a general situation, galaxies 
are described by density profiles which do not allow for analytical calculation 
of the relation r = r{M). Therefore, we must resort to numerical calculations 

dcf 

and take p(M) — p(r(M)), where r(M) is the inverse function to M(r) and 
is found numerically. 

In order to find the areal radius R(t ,M), the constants M, and Ri of 
Eq. (10) have to be appropriately chosen. We may put Ri = 0, thus assuming 
R(to, Mi) — 0. Further we specify Mj to be Ms - the mass already swallowed 
by the singularity. On a spacetime diagram in the coordinates (M, t) that mass 
is represented by the value of M, for which the line t — to, representing the 
present moment, intersects the crunch time function, that is t = tc(Ms) (see 
Fig. 4). Therefore Eq. (10) becomes 



The value of Ms can be calculated using the properties of the apparent horizon 
(see [5], [8]). Only for p{M) given in a fairly simple form analytical calcula- 
tions are possible. Here, we do not find the value of Ms and make use of the 
properties of the apparent horizon. 

The apparent horizon consists of the events for which 



where Mbh is the mass of the black hole [8,30]. We also have Ms < Mbh, 
that is the mass swallowed up by the singularity must necessarily be smaller 
than the mass of the black hole. Therefore we have 



The first term, from the definition of the apparent horizon, is equal to (2Mbh) 3 
and the equation for calculating the value of areal radius for the galaxy under 
consideration becomes 




(12) 



R(M BH ) = 2M BH , 



(13) 




(14) 




(15) 
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Fig. 1 (a) The time of last scattering of CMB in the FLRW model as a function of Qa for 
different values of i?o ■ For the ylCDM model with WM AP 7-year results the instant of last 
scattering is at t ICC = 494 kyr, which is marked with a cross in the plot, (b) The age of the 
Universe for the same set of Qa and Qq values. The cross at t novl = 13.727 Gyr indicates 
the age of the Universe for the /1CDM model with WMAP 7-year results, (c) a = R/M 1 / 3 
and (d) b = R, t /M 1 / 3 at the time of last scattering of CMB in the FLRW model for the 
same set of (2a and i?o values. For the ylCDM model with WMAP 7-year results a equals 
1.649 ■ 10 5 m- 1 / 3 , while b is equal to 3.483 ■ 1(T 3 m" 1 / 3 . 

The advantage of this equation is that it can be applied to any galaxy with 
known density profile and known mass of the black hole. 

4.3 The parameters at the recombination 

In [5] the parameters at the time of last scattering where chosen to match 
those of the k = FLRW model. In this work we will use general FLRW 
models with k ^ and with nonvanishing cosmological constant. The areal 
radius R(t, M) and the scale factor S{t) of the FLRW model are connected by 

R(t, M) — r(M) S(t), (16) 

with M as the radial coordinate. Using the limiting conditions for the FLRW 
model arising from the LT model, i.e. Eq. (8), we have 

R(t, M) = M- l/3 M 1 / 3 S(t), (17) 
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where the constant M§ determines the particular FLRW model. The value of 
this constant can be found by comparison of the solutions of the LT model 
(Eqs. (4)-(6)) with the solutions of the FLRW model [8], obtaining 

M = ~np(t)S 3 (t), (18) 

where p{t) is the (homogeneous) density. By taking the values of p(t) and S(t) 
at the present moment t = t , that is 

p(t ) = n m p crit = n m s(t ) = i, (19) 

where J? m is the density parameter of matter (with equation of state p = 0) , 
p cr it is the critical density with the current estimate of the Hubble constant 
H = 71± 2.5 km s" 1 Mpc" 1 , based on the latest WMAP results [14], we get 
the value for M as 

M = ^Q m Hl (20) 

The value of the scale factor S(t) at the time of recombination (last scat- 
tering of CMB) is found from 

S(t ,„> . JEsL . (21) 

_L i -i-rcc ~ -^rcc 

where z rcc = 1088.2 ± 1.2 is the redshift of CMB from the WMAP results. The 
derivative of the scale factor is determined by the Friedmann equation 

S, t \ 8ttG k , . N 

f J =—*n-g S+ A, (22) 

with the constant k given by 

k = HlS 2 {t n ){Q m + O a -1) = H 2 (f2 Q - 1), (23) 

where Qjy = (3/8ttG)(A/ p CT - lt ) and f2 = fi m + fij\. 

We would like to find the dependence of our model on the fi m and Qa 
parameters of the FLRW model used to determine the functions at the time 
of last scattering. In order to do that we need to know the age of the Universe 
and the recombination time for FLRW models with different fl m and Qa- 
From the Friedmann equation, using the normalization S(to) = 1, the age of 
the Universe to is determined by 



/ ^ . 



The instant of last scattering of CMB, used as one of the input parameters for 
the model, is dependent on the matter content of the Universe. This instant is 
defined by such density of the radiation that the energy of the photon becomes 
smaller than the ionization energy of the hydrogen atom. As the radiation 
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and matter have different equations of state, the instant of last scattering is 
different for different density parameters. The time of last scattering of CMB 
can be found from the equation similar to Eq. (24), but with the upper limit in 
the integral replaced by S(t ICC ), taken from Eq. (21). Figure 1 shows the age of 
the Universe and the instant of recombination as a function of Qa for different 
values of Qq • This figure also shows the dependence of the two functions used 
as input data for the model a(M) = R(M)/M 1 / ?} and b(M) = R, t (Mj/M 1 / 3 
at the recombination on the value of Qa for the same set of Qq values. The 
two parameters are constant within a FLRW model and their values for the 
ylCDM model, in units described in 4.4 and with WMAP 7-year results, are 
a = 1.649 • 10 5 m- 1 / 3 and b = 3.483 • 1(T 3 m" 1 / 3 . 

4.4 The units 

Throughout this article we have used the geometric units where both c and G 
are equal to unity. In order to avoid large numerical values we used 10 11 M Q 
as the unit of mass. This value corresponds roughly to the mass of the M87 
galaxy. 

5 Application of the model to the M87 galaxy 

We have tested the constructed model on the M87 (NGC 4486) galaxy, which 
hosts a supermassive black hole of mass 3.2 • 10 9 M Q and is about 16 Mpc 
away from the Earth [13]. The galaxy is mainly known for its highly relativis- 
tic jet, but we are only interested in its density profile and the mass of the 
supermassive black hole in the center. 

5.1 The M87 input parameters 

In [5] the authors used the approximated profile of Fabricant et. al. [12], with 
the mass of the black hole Mbh = 0.03 m. As was described in Sec. 4.2, we no 
longer have to use any approximation to the density profile and generally we 
could use the original, not approximated profile, and the value of Mbh for our 
numerical calculations. However, in recent years the standard value of Mbh 
and the shape of the density profile for M87 have changed (see Fig. 2(a)). We 
will use the updated values. 

The M87 galaxy is a radio source, and has been the subject of many mea- 
surements based on different methods. In this work we used the 'isotropic' 
profile of Wu, Tremaine [16] based on globular clusters data, given by 



with the values of the parameters: r = 19 kpc, p = 1.9 • 1O 7 M kpc 3 
and a = 1.8. The mass of the black hole is assumed to be 3.2 • 10 9 M Q . The 




(25) 
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Fig. 2 (a) Comparison of M87 density profiles. All curves start at M = Mbh, which for 
Fabricant et. al. is equal to 0.03 m, and 0.032 m for both Wu, Tremaine profiles. The profile 
used in our calculation is 'Wu, Tremaine isotropic'. Both axes are in geometric units, (b) The 
distance-mass relation for the M87 galaxy with the 'Wu, Tremaine isotropic' density profile. 
The mass equal to the mass of the black hole is enclosed within the radius tbh = 25538.6 m. 
(c) R(M) and (d) a(M) functions for the chosen density profile of Eq. (25). The dotted 
vertical line represents the black hole's surface at Mbh = 0.032 m. 



'anisotropic' profile based on the same data, but modified method, differs only 
in the value of a, which for this profile is equal to 1.6. The three density profiles 
are plotted in Fig. 2. Note the change in the magnitude of the density values. 

The quantities used in our calculations that describe the M87 galaxy, be- 
sides Mbh, are the distance-mass relation r(M), the areal radius R(to,M) 
and a(M) = R/M 1 / z . The distance-mass relation is found numerically by in- 
tegration of the density profile and the areal radius R(to,M) is found using 
Eq. (15). The three functions are plotted in Figs. 2(b)-2(d). 



5.2 The scheme of the method 



The general scheme of our method is as follows. Firstly, for the given velocity or 
density profile at the last scattering and the M87 density profile (together with 
the assumed Mbh), the functions E(M) and t B {M) (and also the resulting 
tc(M)) are found in the part outside the black hole, that is for M > M B h- 
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Secondly, for the chosen type of black hole interior, that is the bang and crunch 
time functions for M < Mbh described by the equations in Table 1, random 
four parameters of those equations are chosen and their values are set using 
uniformly distributed random numbers in the appropriate range, depending 
on the parameter. Then, the values of the other four parameters are found by 
numerically solving the equations (11). This is done using multidimensional 
root finding method based on the Hybrid Method. In this way an LT model 
model near M = is found (necessarily recollapsing) , which will serve as the 
starting point for further actions. With some density profiles (at t = now and 
at last scattering) it may turn out that the model is not recollapsing at present. 
This step is necessary to make sure that this does not happen. 

As stated before, for galaxies that host supermassive black holes, such as 
M87, we are mainly interested in finding such an LT model, that would yield 
a black hole of the age of Gyr. Therefore the age of the black hole is the 
parameter to be optimized. Optimization is done using the Simplex algorithm 
of Nelder and Mead for minimization of multidimensional functions [31] by 
again choosing arbitrary four parameters of bang and crunch time functions 
in the interior of the black hole (found in the previous step) and finding such 
values of those parameters that would give a black hole of maximal age. For 
each run we have also found the LT model of minimal age of the black hole, 
just as a comparison to the maximal one. After this we have performed the 
reconstruction of evolution together with the calculation of other valuable 
data. 



5.3 The ACDM model of initial fluctuation 

The current concordance cosmological model is the homogeneous /1CDM model 
with H = 71 ± 2.5 km s" 1 Mpc" 1 , f2 m = 0.266 and Q A = 0.734, based on 
the latest WMAP results [14]. In our calculation we used the parameters of 
this model, described by the equations in Sec. 4.3, for the initial state. The 
constructed model of a galaxy with a black hole takes as the input data for 
the initial fluctuation at last scattering the profiles of either velocity R, t {M), 
or density p{M). Firstly we will present the results of optimization of black 
hole's age for the M87 galaxy evolving from the velocity profile coinciding with 
that of the ylCDM model. 

5.3.1 The velocity profile 

The quantity that is used in our calculations for this evolution type is b = 
R, t (M^/M 1 / 3 . In an FLRW model it depends only on time, and for Q m — 
0.266, Qa = 0.734 its value at recombination is equal to 3.483 • 10~ 3 m~ 1/3 . 
The homogeneity of the chosen velocity profile does not prevent the creation 
of a galaxy as the corresponding profile of density at last scattering, shown in 
Fig. 3(a), is nonuniform, and an initial fluctuation appears and can gain mass 
through accretion. 
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Fig. 3 (a) The ratio of the density at the recombination p rC c(M) corresponding to the 
homogeneous velocity profile of the standard ylCDM model to the ylCDM density in the 
same model. The nonuniformness of this profile is responsible for the creation of a galaxy 
with a central black hole, (b)-(c) The main panel on both plots shows the ig(M) and 
E(M) functions for the maximal and minimal age of the black hole. The maximal age is 
12.719 Gyr, the minimal one is 0.428 Gyr. The functions are calculated for an evolution 
from a homogeneous velocity profile at the recombination consistent with the /1CDM model 
with Q m = 0.266, Qa = 0.734. The bang and crunch time functions inside the black hole 
are parameterized by the hyperbolic function of Table 1. For M > Mbh the functions are 
identical, and dependent only on the input profiles. The difference is inside the black hole 
(M < A/bh) and is exhibited in the inset. The ratio of the two functions is plotted in the 
bottom panel, (d) The main panel shows the tc(M) function for the LT model of maximal 
and minimal age of the black hole. The optimization parameters are the same as described 
in (b) and (c). Note the logarithmic scale for the x axis. The y axis is logarithmic up to the 
value of tc(M^n) = 8.7945 • 10 11 m R3 tg, above that value the scale is linear. This is done 
in order to make the plot readable, as there is a significant difference in the black hole's 
creation time, that is the value of *ah+ a t the minimum, which for the maximal model is 
equal to 6.458 • 10 m. The ratio of the two functions is plotted in the bottom panel. The 
dotted vertical line in all plots represents the boundary of the black hole at Mbh = 0.032 
m. 



Figures 3(b), 3(c) and 3(d) show the main results of optimalization of the 
age of black hole for the evolution from flat velocity profile at recombination. 
The bang and crunch time functions were parameterized by the hyperbolic 
functions of Table 1. The black hole of maximal age is formed at 6.458- 10 10 m, 
what means that the age of such a black hole is 12.719 Gyr. The corresponding 
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quantity 


min 


max 


BH age [Gyr] 


0.428888 


12.719398 


(tAH+)min [10 10 mj 


85.197 


6.458 


61 [10 5 ] 


1.156 


1.160 


f>2 


374.749 


-8.460 ■ 10~ 6 


£>3 


-0.432 


64.205 


64 


0.629 


10.819 


ci [10 11 ] 


4.121 


8.855 


C2 


3.628 • 10 11 


-3.858 ■ 10- 19 


C3 


3.163 


-153.308 


C4 


0.641 


70.526 



Table 2 The parameters of the black hole's bang and crunch time functions for the models 
of minimal and maximal age of the black hole. The parameters labeled as bi describe the 
bang time function, Ci describe the crunch time function. The initial conditions correspond 
to those of Fig. 3. 



values for the minimal age are 8.5197 • 10 11 m and 0.428 Gyr. The last number 
is roughly equal to the value obtained in [5]. Other parameters of the model 
are summarized in Table 2. 

As we can see in Fig. 3(d) the value of the crunch time function at Mbh 
is extremely close to to- However, in the model the crunch time is equal to 
to not for Mbh, but for Ms, the mass swallowed by the singularity, which is 
necessarily smaller than the mass of the black hole. On the other hand Mbh 
is that value for which the future apparent horizon crosses the present time. 
Figure 4 shows this situation in an illustrative way, with the bang and crunch 
time functions chosen so as to make the plot readable. 

Such behaviour is a consequence of the very small difference between the 
future apparent horizon and the crunch time function, which causes Ms and 
Mbh to almost coincide. The value of M5 cannot be found analytically for a 
general observational density profile of the galaxy. In our calculations we have 
used the method described in Sec. 4.2 that does not require the knowledge of 
the value of Ms- 

The maximal age of the black hole obtained here is in good agreement with 
the estimates based on other models. Figure 5(a) shows the time and radial 
dependence of the areal radius R{t, M). As follows from the negative value of 
E, the whole region of spacetime under consideration is recollapsing by now. 
The flat part of the graph, for which t > (tAH+)min, is the region of spacetime 
where the singularity already exists. 

In Fig. 5(b) the constant M slices of R(t, M) for characteristic M values 
are shown. The first nonvanishing M = const slice is for a small M, chosen 
arbitrarily to be 10~ 30 . The whole region under consideration is recollapsing 
by now, but only the region M ^ Mbh is within the black hole. Characteristic 
t = const slices are plotted in Fig. 5(c). The upper plot shows the t = const 
curve for the instant of recombination, the middle one shows the curves for the 
time of the creation of the black hole and for the time it was half its present 
age, the bottom one shows the curve for the present moment. All the curves 
for t values ranging from the recombination time t ICC to the moment of the 
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= 0) 
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time ^--^ 





mass 
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Fig. 4 Evolution leading to a black hole in the E < LT model. This a copy of Fig. 1 from 
[5], where the bang and crunch time functions were chosen so as to make the plot readable 
and illustrative. Intersections of the line t = now with the lines representing the Big Crunch 
and the future apparent horizon AH+ determine the masses Ms and Mbh, respectively. 
AH- stands for past apparent horizon and ths for time, when the future apparent horizon 
first appears. 



black hole creation (tAH+)min start at the point M — 0, R — 0. All the curves 
for larger t values start at M ^ 0, R — 0, because of the formation and growth 
of the black hole. 



5.3.2 The interior parametrization of the black hole 

All the results presented so far were found for the bang and crunch time 
functions of the black hole interior given by the hyperbolic function of Table 1. 
However, the calculations can be performed for an arbitrary type of the interior 
of the black hole. Figure 6 presents the comparison between the hyperbolic and 
the exponential parametrization. The bang and crunch time functions are very 
similar, their ratio differs from unity by a number of the order of 10~ 8 for the 
bang and 10 -3 for the crunch time function. The maximal age of the black 
hole is almost the same for the two types of parametrization. Therefore in the 
following we have chosen only the hyperbolic parametrization of the interior 
of the black hole functions. 
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R(t,M) 




0.001 0.01 M BH 0.1 1 



M [m] 

(c) 

Fig. 5 (a) The function R(t, M) for the evolution from flat velocity profile at the recom- 
bination. The M axis is in logarithmic units in order to enhance the region near the center 
of symmetry. Note the inscribed crunch time function in the R = plane - the crunch time 
function borders the black hole's interior, (b) Characteristic M = const slices of the areal 
radius. Only those slices for which M < Mbh reach the black hole. The M = 1 slice, corre- 
sponding to the region furthermost from the center of symmetry, is recollapsing, but not yet 
in the final singularity. The M = Mbh slice also does not reach the black hole - its value at 
t = to is 2Mbh- 'tree' stands for the time of recombination, (c) Characteristic t = const slices 
of the areal radius. The upper plot shows the curve for the instant of recombination, the 
middle one shows the curves for the time of the creation of the black hole and for the time 
when it was half its present age, the bottom one shows the curve for the present moment. 
The curves for times after the creation of the black hole do not start at M = as the black 
hole is already created. 



5.3.3 The density profile 



An initial fluctuation at the recombination can also be defined by a density 
profile, that is the function p rec (M). The quantity that is used here in calcu- 
lations is a = i?/M 1 / 3 , where the areal radius at the recombination is found 
using the modified Eq. (15), without the first term in square brackets and 
setting the lower limit of the integral to zero. As in the evolution from a given 
velocity profile at last scattering, we set the density profile to the one corre- 
sponding to the /1CDM model with Q m = 0.266 and Q A = 0.734. Then the 
density profile is flat and the value of a is equal to 1.649 • 10 5 m -1 / 3 . 
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Fig. 6 The difference between the hyperbolic and exponential parametrization of the black 
hole interior. The plots show the ratio of the bang (a) and the crunch (b) time functions 
subtracted from unity in order to make the values readable. The difference between the 
functions can be neglected. The shape of the bang time curve for M > Mbh is due to the 
numerical errors, as this M range corresponds to the exterior of the black hole where the 
type of parametrization of the interior does not affect the values of the functions. 



The corresponding profile of velocity R lt (M) is presented in Fig. 7(a). 
This profile exhibits nonuniformness enabling the creation of a black hole. 
However, the nature of the fluctuation is different than in the case of a flat 
velocity profile and its corresponding density profile in Fig. 3(a). In the model 
the density at the center increases by accretion of background matter, which 
is achieved by a lower expansion rate near the center and a higher rate further 
away. The difference in expansion rates is the underlying reason of accretion 
and the creation of black hole. Therefore, any velocity profile at last scattering 
that could lead to the creation of a galaxy with a black hole must necessarily 
be monotonically increasing as a function of M. Such a behaviour can be seen 
in Fig. 7(a). 

The results of our calculations for the evolution form a flat density profile 
coinciding with the ylCDM profile are presented in Figs. 7(b)-7(d). The E(M) 
functions are similar in both cases and the whole M range is recollapsing which 
is necessary for the creation of a black hole. The main difference compared 
to the previous evolution type is in the bang time function ts(M) which is 
increasing. This means that the conditions of avoidance of shell crossings are 
not fulfilled and a shell crossing is imminent. However, due to the values of 
ts{M) it will occur after the present time and not in the range of application 
of the model. Figure 8(a) shows the points where the shell crossings appear. 
We can see that for the points near the black hole it appears after the present 
moment. The non-occurrence of the shell crossing in the area covered by the 
model can also be seen in Fig. 8(b) where there are no points (besides the black 
hole interior) where R,m = 0. The small difference in the LT model functions 
means that the ages of black holes in both evolution scenarios are the same, 
within the numerical errors. Other parameters of the bang and crunch time 
functions are shown in Table 3. 
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quantity 


min 


max 


BH age [Gyr] 


0.429112 


12.719378 


(tAH+)min [10 10 mj 


8.5195 


6.4583 


61 [10 4 ] 


-2.890 


-2.911 


f>2 


-94.851 


2.141 ■ 10~ 6 


£>3 


-0.432 


64.205 


64 


0.629 


10.819 


ci [10 11 ] 


4.119 


8.855 


C2 


3.630 • 10 11 


-3.858 ■ 10- 19 


C3 


3.163 


-153.308 


C4 


0.641 


70.526 



Table 3 The parameters of the bang and crunch time functions of the interior of the black 
hole for the LT models of minimal and maximal age of the black hole for evolution from a 
flat density profile coinciding with the standard /1CDM model. As before, the parameters 
labeled as bi describe the bang time function, Ci describe the crunch time function. 



5.4 The parameters at the recombination 

The results presented so far were obtained for a fiat velocity and flat density 
profiles at the recombination consistent with the /1CDM model with Q m = 
0.266 and Qa — 0.734. However, the specific initial FLRW model should not 
play a key role in the creation of Gyrs old black holes. Therefore, we have 
performed calculations of the maximal age of the black hole for two types 
of evolution for initial models with different parameters. We set the Hubble 
parameter value at Hq = 71 km s _1 Mpc" 1 and the redshift of the CMB at 
z rcc = 1090, based on WMPA 7-year results, and carried out the calculations 
for varying f2 m and Qa values. 

5.4.1 The age of the black hole 

In Fig. 9 the dependence of the maximal age of the black hole on the matter 
density parameter Q m , for three constant values of Q$ equal to 0.5, 1.0 and 1.5, 
is shown. For both evolution types we can clearly see the two constant values of 
the maximal age of the black hole equal to 3.717 and 12.708 Gyr. For Qo = 0.5 
only the first value is obtained. The shapes of the curves for the same value of 
Q an d different evolution types are very similar. The differences are caused by 
numerical errors and are remnants of the method of choosing initial conditions 
for the optimization method involving random numbers. These differences are 
presented in the insets in the Qq = 1.5 plots. 

5.4.2 The LT functions 

Figure 10 shows the E(M) and tc{M) functions for the same three values of 
Qo & s i n Fig. 9 for the evolution form flat velocity at the recombination. For 
the E{M) function all the curves are similar, without any hints of the two 
values of the age of black holes. The plots for the crunch time function reveal 
the existence of the two values, as for all values of Q we can indicate the 
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M [m] M [m] 

(c) (d) 



Fig. 7 (a) The ratio of the velocity at last scattering corresponding to a fiat density profile, 
to the ylCDM velocity. A monotonic increase of this function enables the creation of a galaxy 
with a black hole, (b) The bang time function for the evolution from the flat density profile 
for the maximal and minimal age of the black hole models. The main panel shows the 
functions for the whole M range, the inset shows the functions near the center of symmetry 
at M = and the bottom one shows the ratio of the maximal model to the minimal. Note the 
difference from the bang time function for the evolution from flat velocity profile in Fig. 3(b). 
This function being increasing indicates that there is a shell crossing in the model, but it 
appears outside the range of the application of the model, (c) The corresponding E(M) 
function in the same range as for the velocity evolution in Fig. 3(c). The difference between 
this and the fiat velocity evolution is very small, (d) fn the main panel: the crunch time 
function for the maximal an minimal black hole models. The age of the black hole for the 
maximal model is 12.7f9 Gyr, and 0.429 Gyr for the minimal one, the same values (within 
numerical errors) as for the evolution from a fiat velocity profile. Note the logarithmic x 
axis. 



two distinctive families of curves. Because of the fact that the crunch time 
function and the future apparent horizon t\H+ almost coincide in all the cases 
and taking into account that the value of the future apparent horizon at its 
minimum is the creation time of the black hole, we can say that the crunch 
time curves starting at lower values correspond to the higher value of the 
maximal age of the black hole. The reason why all the crunch time curves 
corresponding to one value of the age of the black hole, e.g. 12.708 Gyr, do 
not start at the same point is that with varying £2 m the age of the Universe 
changes. 
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Fig. 8 (a) The coordinates of the points where the shell crossing singularity occurs. The 
shell crossing appears after the present moment, i.e. not in the range of application of the 
model, (b) Evolution of the areal radius R(t, M) for the evolution from a flat density profile 
at the recombination. There are no shell crossings as there are no points (besides the black 
hole interior) where R,m = 0. 



6 Summary and conclusions 

We have demonstrated that within the inhomogeneous LT model a perturba- 
tion of the velocity or density profile at the recombination can evolve into a 
spherically symmetric galaxy-like object with a black hole at the center of the 
age of Gyr. This work is a refinement to the model presented in [5], which 
allows for such an early creation of black holes, together with the usage of an 
arbitrary density profile describing the present day galaxy and an arbitrary 
FLRW model as the background matter reservoir. The evolution of such a 
perturbation leading to the present mass distribution of a galaxy progresses 
without any shell crossing singularities. In the model, the mass at the center 
increases in consequence of different expansion rates in the central region and 
farther away. The crucial point for the time of the black hole creation is the 
interior geometry of the black hole, that is the bang time function ts and the 
crunch time function tc(M) in the region M < Mbh- By an adequate adjust- 
ment of those arbitrary functions, we were able to find an LT model yielding 
the desired age of Gyrs. It was found that the maximal age of the black hole 
as a function of the density parameters Q m and Qa takes a constant value 
irrespective of the Q$ value. 

The main target for the model was the proper simulation of the black 
hole at the center of M87 galaxy. This supermassive black hole is one of the 
biggest supermassive black holes found. We have used a new density profile 
describing the mass distribution in this galaxy together with an updated value 
of the mass of the black hole based on [13]. This required a slight improvement 
of the basic model allowing to use any density profile, not only the one for 
which analytical calculations are possible. We have successfully obtained an 
LT model characterizing the evolution leading to this galaxy together with 
the black hole created about 12.7 Gyrs ago. This number was independent on 
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Fig. 9 The maximal age of the black hole for three different values of Qq for a flat initial 
velocity profile (left) and a flat initial density profile (right) as a function of Q m ■ The interior 
of the black hole is parameterized using hyperbolic functions. The two constant values of the 
maximal age of the black hole are clearly visible. In the Qq = 0.5 case, for the whole range of 
Q m values, only 3.717 Gyrs old black holes are formed. For the other two values of Qq also 
12.708 Gyrs old black holes are created. In the i?o = 1-0 plots the cross indicates the /1CDM 
model. The insets in the i?rj = 1-5 plots show the small changes of the age values, caused 
by numerical approximations and the usage of random numbers for initial conditions. 



the type of the initial perturbation (velocity or density) seeding the galaxy 
formation. By changing the parameters of the initial FLRW model we found 
that the maximal age of the black hole for this particular galaxy can take only 
two values: 3.717 and 12.708. A set of thorough runs for three values of Q§ 
equal to 0.5, 1.0 and 1.5, determining the three families of FLRW universes: 
open, flat and closed respectively, showed that the LT model functions behave 
properly in all cases. 
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M BH 0.1 0.2 0.3 4 0.5 0.6 0.7 0.8 0.9 1 0.001 0.01 M BH 01 1 

M [m] M [m] 

Fig. 10 The plot shows the E(M) (left) and tc(M) (right) functions for the same three 
values of fio as in Fig. 9. All the plots of E(M) are similar and do not reveal the existence 
of the two values of the age of the black hole. The crunch time functions show this fact by 
the occurrence of the two distinctive families of curves. The curves starting at lower values 
correspond to the higher value of the maximal age of the black hole. Note the logarithmic 
scale on the x axis. The 'flatness' of the curves for M > A/bh is the consequence of the 
logarithmic x axis - in reality the curves are not flat, but change very slowly compared to 
the part inside the black hole. 



However, as found by Gebhardt and Thomas [32], the mass of the black 
hole in M87 can be even 6.4 ± 0.5 • 10 9 M Q , that is twice the mass used 
in our calculations. This can affect the results presented here and we plan to 
perform the simulations with this new value of -Mbh in near future. 

The main disadvantage of the LT model in simulating the evolution and 
creation of galaxies is the absence of rotation. Rotation is a key factor in 
galaxy formation, as it slows down the collapse and produces accretion disks 



26 



Przemyslaw Jacewicz, Andrzej Krasinski 



around black holes. Due to this it plays a much bigger role in the later stages 
of evolution, when the collapse of the central body has started. 
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